library(magrittr)
library(tidyverse)
library(Seurat)
library(readxl)
library(cowplot)
library(colorblindr)
library(viridis)
library(progeny)
library(destiny)

coi <- params$cell_type_super
cell_sort <- params$cell_sort
cell_type_major <- params$cell_type_major
louvain_resolution <- params$louvain_resolution
louvain_cluster <- params$louvain_cluster

1 Cluster markers

1.1 Major T.super markers for cell assign

### load all data ---------------------------------
source("_src/global_vars.R")

# seu_obj <- read_rds(paste0("/work/shah/isabl_data_lake/analyses/16/52/1652/celltypes/", coi, "_processed.rds"))
seu_obj <- read_rds(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/outs_pre/", coi, "_seurat_", louvain_resolution, ".rds"))
# seu_obj <- read_rds(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_highqc.rds"))

myfeatures <- c("umapharmony_1", "umapharmony_2", "sample", louvain_cluster, "doublet", "nCount_RNA", "nFeature_RNA", "percent.mt", "doublet_score")

plot_data_wrapper <- function(cluster_res) {
  cluster_res <- enquo(cluster_res)
  as_tibble(FetchData(seu_obj, myfeatures)) %>% 
    left_join(meta_tbl, by = "sample") %>% 
    rename(cluster = !!cluster_res) %>% 
    mutate(cluster = as.character(cluster),
           tumor_supersite = ordered(tumor_supersite, levels = names(clrs$tumor_supersite)))
}

plot_data <- plot_data_wrapper(louvain_cluster)

1.2 Subtype currated markers

helper_f <- function(x) ifelse(is.na(x), "", x)

markers_v6_super[[coi]] %>% 
  group_by(subtype) %>% 
  mutate(rank = row_number(gene)) %>% 
  spread(subtype, gene) %>% 
  mutate_all(.funs = helper_f) %>% 
  formattable::formattable()
rank CD4.T.CXCL13 CD4.T.naive CD4.T.reg CD4.Th17 CD8.T.CXCL13 CD8.T.cytotoxic CD8.T.ISG Cycling.T.NK dissociated NK.CD56 NK.cytotoxic
1 CD4 CCR7 CD4 CCR6 CCL3 CCL4 IFI6 ASPM BTG1 AREG FCGR3A
2 CD40LG CD4 FOXP3 IL4I1 CCL4L2 CCL5 IFIT1 CENPF DNAJB1 FCER1G FGFBP2
3 CTLA4 CD40LG IL2RA IL7R CD8A CD8A IFIT2 HIST1H4C DUSP1 GNLY KLRD1
4 CXCL13 IL7R TNFRSF4 KLRB1 CRTAM CD8B IFIT3 HMGB2 EGR1 KLRC1 KLRF1
5 FKBP5 KLF2 TRAC LST1 CXCL13 CRTAM ISG15 MKI67 FOS KRT81 PRF1
6 IL6ST LTB LTB FABP5 CST7 MX1 STMN1 FOSB TRDC SPON2
7 ITM2A TCF7 TNFSF13B GZMB DTHD1 MX2 TOP2A HSPA1A TYROBP
8 MAF TPT1 HAVCR2 GZMA RSAD2 TUBA1B HSPA1B XCL1
9 NMB IFNG GZMH TUBB HSPA6 XCL2
10 NR3C1 LAG3 GZMK TYMS JUN
11 PDCD1 MIR155HG GZMM JUNB
12 TNFRSF4 PHLDA1 HLA-DPB1 KLF2
13 TOX2 PTMS ITM2C MT1E
14 TSHZ2 RBPJ KLRG1 MT1X
15 TNFRSF9 TRGC2

1.3 Subtype cluster markers

# marker_tbl <- read_tsv(paste0("/work/shah/isabl_data_lake/analyses/16/52/1652/celltypes/", coi, "_markers.tsv")) %>% 
#   filter(resolution == louvain_resolution)
marker_tbl <- read_tsv(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/outs_pre/", coi, "_markers_", louvain_resolution, ".tsv"))
# marker_tbl <- read_tsv(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_highqc_markers_02.tsv"))

## Hypergeometric test --------------------------------------

test_set <- marker_tbl %>% 
  group_by(cluster) %>% 
  filter(!str_detect(gene, "^RPS|^RPL")) %>% 
  slice(1:50) %>% 
  mutate(k = length(cluster)) %>% 
  ungroup %>%
  select(cluster, gene, k) %>% 
  mutate(join_helper = 1) %>% 
  group_by(cluster, join_helper, k) %>%
  nest(test_set = gene)

markers_doub_tbl <- markers_v6 %>% 
  enframe("subtype", "gene") %>% 
  filter(!(subtype %in% unique(c(coi, cell_type_major)))) %>% 
  unnest(gene) %>% 
  group_by(gene) %>% 
  filter(length(gene) == 1) %>% 
  mutate(subtype = paste0("doublet.", subtype)) %>% 
  bind_rows(tibble(subtype = "Mito.high", gene = grep("^MT-", rownames(seu_obj), value = T)))

ref_set <- markers_v6_super[[coi]] %>% 
  bind_rows(markers_doub_tbl) %>% 
  group_by(subtype) %>% 
  mutate(m = length(gene),
         n = length(rownames(seu_obj))-m,
         join_helper = 1) %>% 
  group_by(subtype, m, n, join_helper) %>%
  nest(ref_set = gene)

hyper_tbl <- test_set %>% 
  left_join(ref_set, by = "join_helper") %>% 
  group_by(cluster, subtype, m, n, k) %>%
  do(q = length(intersect(unlist(.$ref_set), unlist(.$test_set)))) %>%
  mutate(pval = 1-phyper(q = q, m = m, n = n, k = k)) %>%
  ungroup %>%
  mutate(qval = p.adjust(pval, "BH"),
         sig = qval < 0.01)

# hyper_tbl %>% 
#   group_by(subtype) %>% 
#   filter(any(qval < 0.01)) %>%
#   ggplot(aes(subtype, -log10(qval), fill = sig)) +
#   geom_bar(stat = "identity") +
#   facet_wrap(~cluster) +
#   coord_flip()
  
low_rank <- str_detect(unique(hyper_tbl$subtype), "doublet|dissociated")
subtype_lvl <- c(sort(unique(hyper_tbl$subtype)[!low_rank]), sort(unique(hyper_tbl$subtype)[low_rank]))
  
cluster_label_tbl <- hyper_tbl %>% 
  mutate(subtype = ordered(subtype, levels = subtype_lvl)) %>% 
  arrange(qval, subtype) %>%
  group_by(cluster) %>% 
  slice(1) %>% 
  mutate(subtype = ifelse(sig, as.character(subtype), paste0("unknown_", cluster))) %>% 
  select(cluster, cluster_label = subtype) %>% 
  ungroup %>% 
  mutate(cluster_label = make.unique(cluster_label, sep = "_"))

seu_obj$cluster_label <- unname(deframe(cluster_label_tbl)[as.character(unlist(seu_obj[[paste0("RNA_snn_res.", louvain_resolution)]]))])
plot_data$cluster_label <- seu_obj$cluster_label

cluster_n_tbl <- seu_obj$cluster_label %>% 
  table() %>% 
  enframe("cluster_label", "cluster_n") %>% 
  mutate(cluster_nrel = cluster_n/sum(cluster_n))

marker_sheet <- marker_tbl %>% 
  left_join(cluster_label_tbl, by = "cluster") %>% 
  group_by(cluster_label) %>% 
  filter(!str_detect(gene, "^RPS|^RPL")) %>% 
  slice(1:50) %>% 
  mutate(rank = row_number(-avg_logFC)) %>% 
  select(cluster_label, gene, rank) %>% 
  ungroup %>% 
  mutate(cluster_label = ordered(cluster_label, levels = unique(c(names(clrs$cluster_label[[coi]]), sort(cluster_label))))) %>% 
  spread(cluster_label, gene) %>% 
  mutate_all(.funs = helper_f)

marker_tbl_annotated <- marker_tbl %>% 
  left_join(cluster_label_tbl, by = "cluster") %>% 
  left_join(cluster_n_tbl, by = "cluster_label") %>% 
  select(-cluster, -resolution) %>% 
  mutate(cluster_label = ordered(cluster_label, levels = unique(c(names(clrs$cluster_label[[coi]]), sort(cluster_label))))) %>% 
  arrange(cluster_label, -avg_logFC, p_val_adj)

write_tsv(marker_sheet, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_marker_sheet.tsv"))

write_tsv(marker_sheet, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/supplementary_tables/", coi, "_marker_sheet.tsv"))

write_tsv(marker_tbl_annotated, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_marker_table_annotated.tsv"))

write_tsv(marker_tbl_annotated, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/supplementary_tables/", coi, "_marker_table_annotated.tsv"))

formattable::formattable(marker_sheet)
rank CD4.T.naive CD4.T.CXCL13 CD4.T.reg CD4.Th17 CD8.T.cytotoxic CD8.T.CXCL13 CD8.T.ISG NK.CD56 NK.cytotoxic Cycling.T.NK Cycling.T.NK_1 dissociated dissociated_1 dissociated_2 doublet.Fibroblast doublet.Fibroblast_1 doublet.Monocyte doublet.Plasma.cell
1 IL7R CXCL13 TNFRSF4 KLRB1 GZMK CXCL13 IFIT3 GNLY FGFBP2 STMN1 CENPF CCL4L2 HSPA1A KLF2 DCN IGFBP5 HLA-DRA SOX4
2 CCR7 NMB IL2RA IL4I1 CD8A GZMB ISG15 TYROBP FCGR3A TUBA1B ASPM CCL4 HSPA1B CCR7 IGFBP5 TAGLN CST3 PTCRA
3 KLF2 NR3C1 FOXP3 IL7R CD8B CCL4L2 MX1 AREG SPON2 MKI67 MKI67 IFNG MT1X JUNB RBP1 ADIRF CXCL8 MAL
4 EEF1B2 MAF CTLA4 LTB ITM2C MIR155HG IFIT1 KLRC1 PRF1 HIST1H4C HMGB2 FOS DNAJB1 FOS C7 DCN SPP1 MZB1
5 TPT1 FKBP5 LTB LST1 GZMH TNFRSF9 RSAD2 FCER1G KLRF1 TUBB TOP2A FOSB MT1E DUSP1 MEG3 IGFBP7 S100A9 DNTT
6 EEF1A1 IL6ST RTKN2 TNFSF13B CCL5 HAVCR2 IFIT2 TRDC GNLY TOP2A UBE2C TNF HSPA6 SELL CALD1 SPARCL1 S100A8 TFDP2
7 MAL ITM2A BATF CCR6 TRGC2 RBPJ IFI6 XCL1 KLRD1 TYMS CCNB1 JUN HSP90AA1 IL7R IGFBP4 CALD1 HLA-DQA1 CD1E
8 TCF7 TSHZ2 TNFRSF18 CTSH KLRG1 LAG3 MX2 KLRD1 NKG7 CENPF UBE2S EGR1 HSPH1 AREG RARRES2 MGP CD74 STMN1
9 CD40LG CTLA4 SAT1 AQP3 GZMA IFNG IFI44L KRT81 CX3CR1 HMGB2 PTTG1 CCL3L1 HSPE1 EEF1A1 NR2F2 C11orf96 LYZ ARPP21
10 SELL CD40LG TBC1D4 CCL20 CCL4 CCL3 ISG20 XCL2 GZMB ASPM TPX2 NFKBID HSPD1 CD69 SELENOP MYL9 FTL AC084033.3
11 GPR183 PDCD1 TIGIT NFKBIA CRTAM PTMS HERC5 IGFBP2 PLAC8 NUSAP1 STMN1 CD69 HSPB1 BTG2 MDK TIMP3 APOE CDK6
12 LDHB CD4 GADD45A RORA CST7 CD8A SAMD9L CLIC3 CLIC3 PCLAF KPNA2 NR4A2 HSP90AB1 GPR183 SOX4 MDK HLA-DQB1 MAP1A
13 NOSIP LIMS1 TNFRSF1B TNFRSF25 GZMM CRTAM OAS1 IL2RB PLEK HMGN2 CENPE AC020916.1 DNAJA1 PIK3IP1 ADIRF TPM2 MARCKS AC011893.1
14 SNHG8 TNFRSF4 PMAIP1 CEBPD HLA-DPB1 PHLDA1 TNFSF10 CEBPD TYROBP H2AFZ CKS2 EGR2 JUN CD55 MGP NR2F2 AIF1 GLUL
15 PABPC1 RNF19A UGP2 TNFAIP3 DTHD1 FABP5 STAT1 KRT86 PTGDS PCNA TUBA1B DUSP1 CACYBP DNAJB1 EGR1 FBLN1 C1QB ADA
16 NOP53 CORO1B IKZF2 NCR3 PPP1R14B TIGIT EIF2AK2 TXK EFHD2 HIST1H1B HMMR TNFSF9 HSPA8 FKBP5 STAR IGF1 FTH1 AL357060.1
17 LEF1 RBPJ TNFRSF9 SLC4A10 CD3G KRT86 OAS3 CTSW FCER1G DUT DLGAP5 KLF6 FKBP4 TSC22D3 SFRP4 RAMP1 BASP1 GRASP
18 EIF3E CPM ICOS TMIGD2 THEMIS JAML SAMD9 KLRB1 CST7 CLSPN CCNB2 TNFAIP3 CHORDC1 RACK1 CLU SELENOP C1QA CD1B
19 LTB ZBED2 LINC01943 DPP4 DUSP2 CCL5 XAF1 MATK GZMH SMC4 TUBB4B DUSP2 RGS2 LDHB ADAMTS1 RARRES2 C15orf48 MIR181A1HG
20 RACK1 AC004585.1 IL32 TPT1 CD3D LINC01871 GBP1 CCL3 ADGRG1 ATAD2 NUSAP1 BTG2 FOS SARAF TIMP2 LUM FN1 CCDC26
21 NACA TOX2 SOX4 MYBL1 LYAR CXCR6 EPSTI1 CD7 CCL3 SMC2 BIRC5 IER2 ANXA1 CXCR4 TCEAL4 IGFBP6 MNDA JCHAIN
22 UBA52 DUSP4 ARID5B S100A4 TC2N TNIP3 IFI44 NKG7 HOPX TPX2 HMGN2 PPP1R15A DNAJB4 PLAC8 C1R CARMN G0S2 VIPR2
23 TOMM7 AHI1 CD27 CD40LG EOMES PDCD1 PLSCR1 CD63 IGFBP7 TMPO ARL6IP1 NR4A1 PPP1R15A EEF1B2 WFDC2 IGFBP4 APOC1 ID1
24 SOCS3 ICA1 BIRC3 SPOCK2 CXCR6 HLA-DRB1 IFI35 HOPX ZEB2 HELLS CDC20 NFKBIA ZFAND2A TCF7 FHL2 DST SOD2 SOCS2
25 JUNB ARID5B LAYN LINC01871 HLA-DPA1 GAPDH USP18 TMIGD2 PRSS23 UBE2C CDKN3 JUND FOSB TPT1 IFITM3 CAV1 NPC2 RCAN1
26 SERINC5 CD84 CORO1B ERN1 HLA-DRB1 FAM3C OASL CMC1 AKR1C3 NASP SMC4 GADD45B DNAJA4 FOSB PEG3 ID4 LST1 CLDN5
27 TMEM123 CCDC50 TYMP JAML SLF1 CTLA4 MT2A TNFRSF18 CD247 DEK CKS1B TSC22D3 DUSP1 NOSIP C1S NUPR1 GSN CASC15
28 EEF2 IGFL2 DUSP4 FKBP11 APOBEC3G SPRY1 HELZ2 KLRC2 MYBL1 RRM2 KIF20B TAGAP TSPYL2 SC5D LUM CSRP2 MS4A6A PFKFB2
29 FXYD5 RGS1 CD4 IFNGR1 PPP2R5C CCND2 TRIM22 SRGAP3 AREG CKS1B GTSE1 ZFP36 SERPINH1 NACA SERPINF1 CDKN1C IL1B GALNT2
30 TSHZ2 BATF ENTPD1 MGAT4A KIAA1551 CD63 CMPK2 GSTP1 C1orf21 KNL1 TUBA1C ZFP36L1 UBC AP3M2 CEBPD PGR PSAP HES4
31 TRABD2A SRGN CTSC S100A6 F2R CD8B LY6E LAT2 S1PR5 HMGB1 TUBB RGCC UBB BTG1 AKAP12 COL6A1 GRN MARCKSL1
32 ANK3 CH25H MIR4435-2HG ELK3 SLAMF7 GZMH PARP14 GZMB CTSW MCM7 SGO2 IL7R AHSA1 NOP53 TIMP1 COL6A2 CD83 TP53INP1
33 SARAF SPOCK2 LINC02099 EEF1A1 CXCR4 ENTPD1 NT5C3A LINC00996 ABHD17A FABP5 H2AFZ CRTAM NEU1 ZBTB16 CST3 EMX2 CTSH APBA2
34 AQP3 ZNRF1 MAGEH1 KIT SH2D1A SRGAP3 DDX58 PRF1 KLF2 UBE2S JPT1 KDM6B DNAJB6 ERAP2 NR2F1 PALLD GLUL TSHR
35 RIPOR2 CHN1 SPOCK2 LTC4S FAM102A VCAM1 IRF7 KLRF1 PTPN12 GAPDH CEP55 NR4A3 GADD45B CCND3 DLK1 PPP1R14A MEF2C NREP
36 AP3M2 TNFRSF25 PHACTR2 RUNX2 CD52 ID2 RNF213 NCAM1 TTC38 CXCL13 NUF2 ATF3 KLF6 EEF1D SERPING1 RBP1 CYBB MME
37 ZFAS1 CD200 CARD16 RORC YBX3 HLA-DRA DDX60L CXXC5 KLRB1 RANBP1 KIF14 DDX3X BTG2 PPP1R15A BEX3 C7 CTSB AC002454.1
38 LINC02273 METTL8 S100A4 ZBTB16 STK17A ITGAE DDX60 SH2D1B CCL4 H2AFX KNL1 DUSP6 JUNB PLK3 C11orf96 MFGE8 CSF3R CHI3L2
39 EIF4B RILPL2 STAM FAM241A CCR5 DUSP4 SAT1 MCTP2 PTGDR MCM3 CENPA GZMK TNF ZFP36L2 FILIP1L KANK2 CD14 SMIM3
40 TOB1 TNFRSF18 GLRX PDE4D GPR174 LYST PPM1K IFITM3 ITGB2 EZH2 CDK1 KLRG1 ERN1 EEF2 RARRES1 NR2F1 C1QC SSBP2
41 SESN3 SLA SPATS2L IL23R COTL1 TNFSF4 PNPT1 ZNF683 XBP1 TUBB4B HMGB3 JUNB JUND HNRNPA1 COL1A2 MFAP4 SGK1 UHRF1
42 NSA2 SMCO4 AC005224.3 B3GALT2 CCL4L2 NDFIP2 PARP9 ITGA1 CEP78 H2AFV CDCA8 RASGEF1B ATF3 VSIR NUPR1 COL1A1 SPI1 LRRC28
43 PASK BTLA MAF CERK CD3E AKAP5 IFIH1 IFITM2 ARL4C CDK1 PLK1 ANXA1 DEDD2 PASK TCEAL9 PBX1 FCGRT BCL11A
44 TNFRSF25 NAP1L4 PBXIP1 TLE1 TUBA4A GOLIM4 SP110 CCL5 BIN2 HIST1H1D AURKA MCL1 CD69 EIF3H MARCKSL1 GSN EGR1 SCAI
45 FAU FYB1 F5 EEF1B2 PECAM1 CD27 OAS2 ITGAX LITAF HNRNPAB TROAP CXCR4 CLK1 TXNIP SLC40A1 PDGFRB FCGR2A ATP6AP1L
46 EEF1D MIR155HG SLAMF1 CFH ARAP2 SNAP47 C19orf66 CD38 TRDC CENPE H2AFV IER5 IER5L LDLRAP1 CFH SERPINF1 PLAUR RUFY3
47 LDLRAP1 PTPN13 BTG3 PERP ITGA1 RGS1 STAT2 SAMD3 TXK TK1 KIF2C MYADM H3F3B CMTM8 APOE SFRP1 CPVL CD79A
48 CTSL SESN3 TRAC PLAT JAML HLA-DPA1 LAG3 SLC16A3 MYOM2 ZWINT MAD2L1 CD8A NR4A1 SCML1 GNG11 LHFPL6 MS4A7 GNA15
49 ITGA6 BIRC3 IL1R1 KIF5C CD84 LINC02446 LAP3 CAPN12 GZMM BIRC5 NUCKS1 PTGER4 CXCR4 LINC00402 CDKN1C PLAC9 ALDH2 HHIP-AS1
50 PFDN5 TP53INP1 DNPH1 PLCB1 AOAH SAMSN1 APOL6 CD247 CD300A PTTG1 RAD21 ZFP36L2 EIF4A2 RIPK2 NBL1 SERPING1 SERPINA1 GSTM3

2 Clusters

2.1 sizes

enframe(sort(table(seu_obj$cluster_label))) %>% 
  mutate(name = ordered(name, levels = rev(name))) %>% 
  ggplot() +
  geom_bar(aes(name, value), stat = "identity") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  labs(y = c("#cells"), x = "cluster")

2.2 UMAP

alpha_lvl <- ifelse(nrow(plot_data) < 20000, 0.2, 0.1)
pt_size <- ifelse(nrow(plot_data) < 20000, 0.2, 0.05)

common_layers_disc <- list(  
  geom_point(size = pt_size, alpha = alpha_lvl),
  NoAxes(),
  guides(color = guide_legend(override.aes = list(size = 2, alpha = 1))),
  labs(color = "")
)

common_layers_cont <- list(  
  geom_point(size = pt_size, alpha = alpha_lvl),
  NoAxes(),
  scale_color_gradientn(colors = viridis(9)),
  guides(color = guide_colorbar())
)

ggplot(plot_data, aes(umapharmony_1, umapharmony_2, color = cluster_label)) + 
  common_layers_disc +
  #facet_wrap(~therapy) +
  ggtitle("Sub cluster")

3 Filtering out doublet clusters

3.1 UMAP

my_subtypes <- names(clrs$cluster_label[[coi]])
my_subtypes <- c(my_subtypes, unlist(lapply(paste0("_", 1:3), function(x) paste0(my_subtypes, x)))) %>% .[!str_detect(., "doublet|dissociated")]
my_subtypes <- my_subtypes[my_subtypes %in% unique(seu_obj$cluster_label)]
my_subtypes <- my_subtypes[my_subtypes %in% names(clrs$cluster_label[[coi]])]

cells_to_keep <- colnames(seu_obj)[seu_obj$cluster_label %in% my_subtypes]
# seu_obj_sub <- subset(seu_obj, cells = cells_to_keep)
# seu_obj_sub <- RunUMAP(seu_obj_sub, dims = 1:50, reduction = "harmony", reduction.name = "umapharmony", reduction.key = "umapharmony_")
# seu_obj_sub$cluster_label <- seu_obj$cluster_label[colnames(seu_obj) %in% colnames(seu_obj_sub)]
# write_rds(seu_obj_sub, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered.rds"))
seu_obj_sub <- read_rds(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered.rds"))

plot_data_sub <- as_tibble(FetchData(seu_obj_sub, c(myfeatures, "cluster_label"))) %>% 
  left_join(meta_tbl, by = "sample") %>% 
  mutate(patient_id = str_remove_all(patient_id, "SPECTRUM-OV-"),
         tumor_supersite = ordered(tumor_supersite, levels = names(clrs$tumor_supersite))) %>% 
  mutate(cell_id = colnames(seu_obj_sub),
         cluster_label = ordered(cluster_label, levels = my_subtypes),
         )
  
if (cell_sort == "CD45+") {
plot_data_sub <- filter(plot_data_sub, sort_parameters != "singlet, live, CD45-", !is.na(tumor_supersite))
}

if (cell_sort == "CD45-") {
plot_data_sub <- filter(plot_data_sub, sort_parameters != "singlet, live, CD45+", !is.na(tumor_supersite))
}

ggplot(plot_data_sub, aes(umapharmony_1, umapharmony_2, color = cluster_label)) + 
  common_layers_disc +
  ggtitle("Sub cluster") +
  #facet_wrap(~cluster_label) +
  scale_color_manual(values = clrs$cluster_label[[coi]])

ggplot(plot_data_sub, aes(umapharmony_1, umapharmony_2, color = patient_id_short)) + 
  common_layers_disc +
  ggtitle("Patient") +
  # facet_wrap(~therapy) +
  scale_color_manual(values = clrs$patient_id_short)

ggplot(plot_data_sub, aes(umapharmony_1, umapharmony_2, color = tumor_supersite)) + 
  # geom_point(aes(umapharmony_1, umapharmony_2), 
  #            color = "grey90", size = 0.01, 
  #            data = select(plot_data_sub, -tumor_supersite)) +
  common_layers_disc +
  ggtitle("Site") +
  # facet_wrap(~therapy) +
  scale_color_manual(values = clrs$tumor_supersite)

write_tsv(select(plot_data_sub, cell_id, everything(), -umapharmony_1, -umapharmony_2, -contains("RNA_")), paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_embedding.tsv"))

4 sub clustering

big_clusters <- list(
  CD4.T.CXCL13 = "CD4.T.CXCL13", 
  CD4.T.naive = "CD4.T.naive", 
  CD4.T.reg = "CD4.T.reg", 
  CD4.Th17 = "CD4.Th17", 
  CD8.T.CXCL13 = "CD8.T.CXCL13", 
  CD8.T.cytotoxic = "CD8.T.cytotoxic", 
  CD8.T.ISG = "CD8.T.ISG", 
  Cycling.T.NK = "Cycling.T.NK", 
  NK.CD56 = "NK.CD56", 
  NK.cytotoxic = "NK.cytotoxic"
)

cells_to_keep_list <- lapply(big_clusters, function(x) 
colnames(seu_obj_sub)[seu_obj_sub$cluster_label %in% x])

# seu_list <- lapply(cells_to_keep_list, function(x) subset(seu_obj_sub, cells = x))
# 
# preprocess_wrapper <- . %>%
#   FindNeighbors(reduction = "harmony", dims = 1:50) %>%
#   FindClusters(res = 0.1) %>%
#   FindClusters(res = 0.2) %>%
#   FindClusters(res = 0.3)
# 
# seu_list <- lapply(seu_list, preprocess_wrapper)
# write_rds(seu_list, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered_sub_clusters.rds"))

seu_list <- read_rds(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered_sub_clusters.rds"))

get_cluster_vector <- function(cluster_label, cluster_res = RNA_snn_res.0.3) {
  cluster_res <- enquo(cluster_res)
  cbind(cell_id = colnames(seu_list[[cluster_label]]), FetchData(seu_list[[cluster_label]], c(as_label(cluster_res)))) %>%
    as_tibble %>%
    mutate(cluster_extended = paste0(cluster_label, "_", !!cluster_res)) %>%
    select(cell_id, cluster_extended) %>% 
    deframe
}

cluster_vector_list <- lapply(big_clusters, get_cluster_vector)

cluster_vector <- unlist(cluster_vector_list, use.names = F) %>% setNames(lapply(cluster_vector_list, names) %>% unlist(use.names = F))

seu_obj$cluster_extended <- cluster_vector[seu_obj$cell_id] %>%
  setNames(seu_obj$cell_id)

seu_obj_sub$cluster_extended <- cluster_vector[seu_obj_sub$cell_id] %>%
  setNames(seu_obj_sub$cell_id)

cluster_extended_uniq <- as.character(na.omit(unique(seu_obj_sub$cluster_extended)))

Idents(seu_obj_sub) <- seu_obj_sub$cluster_extended
Idents(seu_obj) <- seu_obj$cluster_extended

# marker_tbl_extended <- lapply(
#   cluster_extended_uniq,
#   function(x) {
#     marker_tbl <- FindMarkers(seu_obj_sub, ident.1 = x) %>%
#       as_tibble(rownames = "gene") %>%
#       separate(gene, into = c("gene", "drop"), sep = "\\.\\.\\.") %>%
#       select(-drop) %>% 
#       mutate(cluster_extended = x)
#     write_tsv(marker_tbl, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_marker_table_subcluster_", x, ".tsv"))
#     return(marker_tbl)
# }) %>% bind_rows
#  
# write_tsv(marker_tbl_extended, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_marker_table_annotated_extended.tsv"))

marker_tbl_extended <- read_tsv(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_marker_table_annotated_extended.tsv"))
 
cluster_label_extended <- c(
  CD4.T.naive_0 = "CD4.T.naive.1", 
  CD4.T.naive_1 = "CD4.T.naive.2", 
  CD4.T.CXCL13_0 = "CD4.T.CXCL13.early", 
  CD4.T.CXCL13_1 = "CD4.T.CXCL13.late.1", 
  CD4.T.CXCL13_2 = "CD4.T.CXCL13.late.2", 
  CD4.T.reg_0 = "CD4.T.reg.2", 
  CD4.T.reg_1 = "CD4.T.reg.1", 
  CD4.T.reg_2 = "CD4.T.reg.ISG", 
  CD4.T.reg_3 = "CD4.T.reg.3", 
  CD4.Th17_0 = "CD4.Th17.1", 
  CD4.Th17_1 = "CD4.Th17.2", 
  CD8.T.cytotoxic_0 = "CD8.T.naive.2", 
  CD8.T.cytotoxic_1 = "CD8.T.cytotoxic", 
  CD8.T.cytotoxic_2 = "dissociated_3", 
  CD8.T.cytotoxic_3 = "CD8.T.naive.1", 
  CD8.T.cytotoxic_4 = "gd.T.cell", 
  CD8.T.cytotoxic_5 = "doublet.Plasma.cell_1", 
  CD8.T.CXCL13_0 = "CD8.T.CXCL13.early", 
  CD8.T.CXCL13_1 = "CD8.T.CXCL13.late", 
  CD8.T.CXCL13_2 = "NK.reg.CRTAM", 
  CD8.T.ISG_0 = "CD8.T.ISG", 
  CD8.T.ISG_1 = "CD4.T.ISG", 
  CD8.T.ISG_2 = "dissociated_4", 
  NK.CD56_0 = "NK.reg.IGFBP2",
  NK.CD56_1 = "NK.reg.CD56", 
  NK.CD56_2 = "NK.reg.ISG", 
  NK.CD56_3 = "NK.reg.KRT81.KRT86", 
  NK.CD56_4 = "NK.reg.CCL3",
  NK.cytotoxic_0 = "NK.cytotoxic.SPON2.2", 
  NK.cytotoxic_1 = "NK.cytotoxic.GZMH", 
  NK.cytotoxic_2 = "NK.cytotoxic.SPON2.1",
  Cycling.T.NK_0 = "Cycling.CD8.T.1",
  Cycling.T.NK_1 = "Cycling.CD8.T.2", 
  Cycling.T.NK_2 = "Cycling.CD8.T.3", 
  Cycling.T.NK_3 = "Cycling.CD8.T.4", 
  Cycling.T.NK_4 = "Cycling.NK.3", 
  Cycling.T.NK_5 = "Cycling.CD4.T", 
  Cycling.T.NK_6 = "Cycling.NK.2", 
  Cycling.T.NK_7 = "Cycling.NK.1"
)

names(seu_obj$cluster_label) <- colnames(seu_obj)
seu_obj$cluster_label[is.na(seu_obj$cluster_label)] <- "NA"
Idents(seu_obj) <- seu_obj$cluster_label
seu_obj$cluster_label_sub <- cluster_label_extended[seu_obj$cluster_extended]
names(seu_obj$cluster_label_sub) <- colnames(seu_obj)
seu_obj$cluster_label_sub[is.na(seu_obj$cluster_label_sub)] <- "NA"
Idents(seu_obj) <- seu_obj$cluster_label_sub

names(seu_obj_sub$cluster_label) <- colnames(seu_obj_sub)
seu_obj_sub$cluster_label[is.na(seu_obj_sub$cluster_label)] <- "NA"
Idents(seu_obj_sub) <- seu_obj_sub$cluster_label
seu_obj_sub$cluster_label_sub <- cluster_label_extended[seu_obj_sub$cluster_extended]
names(seu_obj_sub$cluster_label_sub) <- colnames(seu_obj_sub)
seu_obj_sub$cluster_label_sub[is.na(seu_obj_sub$cluster_label_sub)] <- "NA"
Idents(seu_obj_sub) <- seu_obj_sub$cluster_label_sub


marker_sheet_extended <- marker_tbl_extended %>%
  mutate(cluster_label_sub = cluster_label_extended[cluster_extended]) %>% 
  group_by(cluster_label_sub) %>%
  mutate(rank = row_number()) %>% 
  slice(1:50) %>% 
  select(rank, gene, cluster_label_sub) %>% 
  arrange(cluster_label_sub) %>% 
  spread(cluster_label_sub, gene)

# marker_sheet_joined <- marker_sheet %>% 
#   select(-CD8.T.ISG) %>% 
#   left_join(marker_sheet_extended, by = "rank") %>%
#   gather(cluster_label, gene, -rank) %>% 
#   mutate(cluster_label = ordered(cluster_label, levels = unique(c(names(clrs$cluster_label[[coi]]), sort(cluster_label))))) %>% 
#   spread(cluster_label, gene)

cluster_n_tbl_full <- seu_obj$cluster_label_sub %>% 
  table() %>% 
  enframe("cluster_label_sub", "cluster_n") %>% 
  mutate(cluster_nrel = cluster_n/sum(cluster_n)) %>% 
  filter(cluster_label_sub != "NA")

marker_tbl_extended_annotated <- marker_tbl_extended %>% 
  mutate(cluster_label_sub = cluster_label_extended[cluster_extended]) %>%
  select(-cluster_extended) %>% 
  left_join(cluster_n_tbl_full, by = "cluster_label_sub")

# marker_tbl_annotated_full <- marker_tbl_annotated %>% 
#   filter(!(cluster_label %in% unlist(big_clusters))) %>% 
#   bind_rows(marker_tbl_extended_annotated) %>% 
#   mutate(cluster_label = ordered(cluster_label, levels = unique(c(names(clrs$cluster_label[[coi]]), sort(cluster_label))))) %>% 
#   arrange(cluster_label, -avg_logFC)


cells_to_keep <- colnames(seu_obj_sub)[!(str_detect(seu_obj_sub$cluster_label_sub, "dissociated|doublet"))]
# seu_obj_sub_sub <- subset(seu_obj_sub, cells = cells_to_keep)
# seu_obj_sub_sub <- RunUMAP(seu_obj_sub_sub, dims = 1:50, reduction = "harmony", reduction.name = "umapharmony", reduction.key = "umapharmony_")
# write_rds(seu_obj_sub_sub, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered_sub.rds"))

seu_obj_sub_sub <- read_rds(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered_sub.rds"))

formattable::formattable(marker_sheet_extended)
rank CD4.T.CXCL13.early CD4.T.CXCL13.late.1 CD4.T.CXCL13.late.2 CD4.T.ISG CD4.T.naive.1 CD4.T.naive.2 CD4.T.reg.1 CD4.T.reg.2 CD4.T.reg.3 CD4.T.reg.ISG CD4.Th17.1 CD4.Th17.2 CD8.T.CXCL13.early CD8.T.CXCL13.late CD8.T.cytotoxic CD8.T.ISG CD8.T.naive.1 CD8.T.naive.2 Cycling.CD4.T Cycling.CD8.T.1 Cycling.CD8.T.2 Cycling.CD8.T.3 Cycling.CD8.T.4 Cycling.NK.1 Cycling.NK.2 Cycling.NK.3 dissociated_3 dissociated_4 doublet.Plasma.cell_1 gd.T.cell NK.cytotoxic.GZMH NK.cytotoxic.SPON2.1 NK.cytotoxic.SPON2.2 NK.reg.CCL3 NK.reg.CD56 NK.reg.CRTAM NK.reg.IGFBP2 NK.reg.ISG NK.reg.KRT81.KRT86
1 CXCL13 CXCL13 CXCL13 IFIT3 CCR7 IL7R RTKN2 TNFRSF4 TNFRSF4 ISG15 KLRB1 LST1 CXCL13 CXCL13 GZMH IFIT3 CD8B GZMK STMN1 STMN1 STMN1 CXCL13 CXCL13 SPON2 CENPF CENPF GZMK XIST IGHG3 TRGC2 FGFBP2 SPON2 FGFBP2 CCL3 GNLY XCL2 TYROBP ISG15 GNLY
2 MAF NMB PDCD1 ISG15 RPS8 CD40LG FOXP3 IL2RA TNFRSF18 SAT1 IL7R IL4I1 CCL4L2 KRT86 CD8A ISG15 IL7R CD8A MKI67 TUBA1B MKI67 STMN1 HLA-DRA FCGR3A ASPM MKI67 MALAT1 MALAT1 IGHG1 XCL1 GZMH FGFBP2 FCGR3A AREG AREG GZMB KLRC1 GNLY KLRC1
3 IL6ST CPM CTLA4 MX2 RPS6 ANXA1 LTB FOXP3 TNFRSF9 HERC5 IL4I1 TNFSF13B HAVCR2 TNFRSF9 CD8B IFIT1 CD8A CD8B TUBA1B HIST1H4C TOP2A TYMS HLA-DRB1 IGFBP7 MKI67 STMN1 PTPRC NEAT1 IGLC3 ZNF683 CX3CR1 FCGR3A SPON2 CCL4L2 KLRC1 CRTAM FCER1G IFIT3 KRT81
4 NR3C1 NR3C1 MAF IFIT2 RPS3A GPR183 UGP2 BATF IL2RA MX1 LTB PCDH9 MIR155HG GZMB CCL5 MX1 YBX3 KLRG1 HMGB2 MKI67 ASPM TUBA1B GZMB KLRF1 TOP2A TOP2A EDF1 MT-ND3 IGLC2 KLRC2 NKG7 CCL3 KLRF1 XCL2 ITGAX XCL1 TRDC FCER1G FCER1G
5 FKBP5 IGFL2 ZBED2 IFI44L SELL KLF2 TBC1D4 TNFRSF18 BATF TNFRSF4 CCR6 SPINK2 PTMS RBPJ HLA-DRB1 IFIT2 ANXA1 ITM2C TOP2A TUBB TYMS PCNA IFNG CX3CR1 KNL1 ASPM RPS29 MT-ND1 IGHG4 LINC02446 TRGC2 PRF1 PRF1 XCL1 FCER1G FABP5 GNLY MX1 AREG
6 CD40LG FKBP5 RBPJ IFIT1 EEF1B2 S100A4 IL2RA LTB CTLA4 SPATS2L TNFSF13B KIT CD8A SRGAP3 GZMA RSAD2 RPS12 CCL4 TUBB TYMS HELLS PCLAF FABP5 AKR1C3 NUSAP1 TUBA1B RPL36AL MT-CO3 IGKC HOPX KLF2 KLRF1 GNLY FOS XCL1 TNFRSF9 AREG IFIT1 TNFRSF18
7 CD4 ITM2A DUSP4 RSAD2 RPS12 FOS PMAIP1 CTLA4 GADD45A IFI6 LST1 SCN1B LAG3 FAM3C GZMK IFI6 TOB1 DUSP2 CENPF TOP2A CLSPN DUT COTL1 GINS2 CENPE HMGB2 RPL13A MT-ND2 JCHAIN XCL2 KLRG1 PTGDS KLRD1 IER2 TYROBP HSP90AB1 IGFBP2 TYROBP KRT86
8 AC004585.1 ICA1 CD4 MX1 RPL34 LTB IKZF2 GADD45A TNFRSF1B IL2RA AQP3 LTC4S GZMB TIGIT HLA-DPB1 ISG20 FOS PPP1R14B HMGN2 HMGB2 PCNA FABP5 GAPDH FGFBP2 TPX2 GNLY ATP5MC2 MT-CYB IGHGP CTSW PLEK CX3CR1 NKG7 CCL4 TRDC PKM KRT81 IFI27 TRDC
9 CORO1B IL6ST CD40LG IL7R RPL32 TPT1 IL32 RTKN2 ICOS CTLA4 CCL20 IL1R1 IFNG SPRY1 PTMS MX2 VIM TRGC2 PTTG1 NUSAP1 UBE2C CLSPN NASP MLC1 DLGAP5 TUBB CALM1 MT-ATP6 IGHG2 CD63 PRF1 GZMB CLIC3 EGR1 XCL2 GAPDH KLRD1 IFITM3 TYROBP
10 ITM2A TSHZ2 MIR155HG IFI6 KLF2 NFKBIA CTLA4 SOX4 PIM3 OAS1 CTSH LINC00299 CCL3 IFNG ZNF683 SAMD9L JUNB CST7 IL2RA CENPF PCLAF HELLS LAG3 S1PR5 GTSE1 NUSAP1 SRP14 MT-CO1 CD8A CD7 ZEB2 IGFBP7 TYROBP NFKBIA KLRD1 NME1 CLIC3 KLRC1 XCL1
11 TSHZ2 TNFRSF4 PRDM1 TNFSF10 RPL9 TIMP1 TIGIT TNFRSF1B FOXP3 IFIT3 RORA GSN CXCR6 MIR155HG HLA-DPA1 IFI44L GZMK DTHD1 PCLAF PCLAF GTSE1 NASP MCM5 TYMS KIF14 UBE2C COMMD6 PTPRC MZB1 TRGC1 PRSS23 PLAC8 GZMB CCL3L1 CMC1 MIR155HG KLRB1 AREG CTSW
12 PDCD1 AHI1 TNFRSF18 STAT1 IL7R S100A10 ARID5B TIGIT IL1R2 IFIT1 TNFRSF25 PLAT HLA-DRB1 CD8A CXCR6 HERC5 CCR7 CCL4L2 NUSAP1 HIST1H1B ORC6 MCM3 PTMS E2F1 HJURP IGFBP2 RPS25 MT-CO2 CCL5 NCR3 KLRD1 GNLY PLAC8 DUSP2 CTSW SRM CEBPD IFI44L ZNF683
13 TOX2 AC004585.1 ITM2A HERC5 RPS13 JUN SAT1 TBC1D4 LTB FOXP3 SLC4A10 CD300LF CCL5 LAG3 HLA-DRA OAS1 ZFP36L2 CRTAM ASPM HMGN2 CDT1 MKI67 MIR155HG MCM7 POLQ TPX2 PFDN5 H3F3B IGHA1 CD52 FCGR3A FCER1G CX3CR1 BTG2 NCAM1 NPW KRT86 RSAD2 XCL2
14 CH25H LIMS1 RGS1 ISG20 RPL21 S100A11 BATF LINC01943 TBC1D4 TYMP DPP4 AFF3 RBPJ PHLDA1 ITM2C MT2A S100A6 GZMM H2AFZ ASPM BRCA1 MCM5 HAVCR2 MYOM2 BUB1B KLRC1 RPS15A CCNI CCL4 TRDC MYBL1 NKG7 FCER1G NR4A1 HOPX KDM6B TXK CD38 HOPX
15 PTPN13 BIRC3 LINC01943 OAS1 TCF7 VIM CORO1B SAT1 PKM LTB MYBL1 DLL1 GZMH ENTPD1 COTL1 STAT1 MT-ND1 CMC1 TYMS H2AFZ ZWINT MCM7 CD8A MCM4 HMMR HMGN2 HINT1 GMFG GZMK KLRD1 ADGRG1 AKR1C3 PTGDS FOSB LINC00996 TSPAN17 XCL1 IFIT2 RGS16
16 CTSW PGM2L1 TNFRSF4 OAS3 MAL CTSL MIR4435-2HG TNFRSF9 BIRC3 RTKN2 S100A4 LIF HLA-DRA HAVCR2 HLA-DQA1 TNFSF10 KLF2 LYAR CEP55 RRM2 ATAD2 MCM4 HELLS CHEK1 KIF11 SMC4 RPS27 COX5B CD8B KIR3DL2 EFHD2 AREG PLEK CD69 KLRF1 ZBED2 IL2RB IFI6 CD7
17 CCL5 CD200 ADAM19 EIF2AK2 RPL39 RPL34 TNFRSF4 UGP2 REL EPSTI1 CD40LG IL23R PHLDA1 CCL3 IFNG OASL LYAR CCL5 DUT SMC4 RAD51AP1 TUBB MCM6 CLSPN CIT FCER1G FTH1 CALM2 RGS1 TRG-AS1 FLNA CCL4 EFHD2 TYROBP IFITM3 GOLIM4 MATK TRDC LGALS1
18 NKG7 SESN3 CCDC50 IFI44 RPL22 EEF1A1 AC005224.3 IKZF2 ENTPD1 STAT1 TPT1 RORC PDCD1 CTLA4 CCL4 PLSCR1 CXCR4 EOMES PCNA PCNA NUSAP1 GAPDH ACTB STMN1 SGO2 TYMS HSPA8 PSMB1 GZMA CCL5 GNLY KLRD1 CST7 NR4A2 TXK PAICS NKG7 TNFSF10 CD63
19 KLRD1 RNF19A TSHZ2 GBP1 RPL7 RPLP0 PBXIP1 PMAIP1 TYMP ICOS RPLP1 TGM2 HLA-DPA1 JAML APOBEC3G GBP1 NELL2 SH2D1A TNFRSF4 DUT DLGAP5 GINS2 MCM4 CHST2 C21orf58 GSTP1 DBI TPM3 HLA-DPA1 TMSB4X CST7 PLEK CCL3 GADD45B IL2RB MRTO4 XCL2 ISG20 KLRC2
20 ST8SIA1 CHN1 CD2 CCR7 RPL5 RPS18 STAM IL32 SAT1 ISG20 EEF1A1 TNFSF11 CD8B CCND2 SH3BGRL3 SAMD9 S100A10 YBX3 TPX2 ATAD2 CDK1 HMGB2 STMN1 CTBP2 KIF23 KRT81 FAU SOD1 CXCR6 LINC01871 HOPX CLIC3 IGFBP7 KLRD1 KLRC2 GEM CD63 GZMB TMIGD2
21 LIMS1 GNG4 ETV7 EPSTI1 RPS18 RPS12 GADD45A ICOS IL1R1 MX2 B3GALT2 ZG16B HLA-DPB1 PTMS CD52 EIF2AK2 MT-ND2 GZMA FOXP3 CLSPN BIRC5 RANBP1 PCNA FCER1G AC007240.1 TRDC PARK7 ATP5MPL IKZF3 GNLY C12orf75 EFHD2 ADGRG1 TRDC MATK EGR2 TMIGD2 SAMD9L SRGAP3
22 SMCO4 MAF MYO7A MT2A RPL37 RPL39 BTG3 LAYN SDC4 TIGIT RORC ENPP1 VCAM1 LINC01871 CD3D LAG3 TNFRSF18 TC2N ENO1 UBE2C TUBA1B HLA-DRA CLSPN PLAC8 KIF15 H2AFZ RPLP2 AES LTB ZFP36L2 GZMB MYOM2 HOPX KLRC1 SELL CD72 CMC1 BST2 MATK
23 CD8A BATF TOX2 SAMD9L RPL13 JUNB BIRC3 LINC02099 LINC01943 RSAD2 LTK PTGER3 LINC01871 GOLIM4 LAG3 IFI35 LMNA CXCR4 BATF HMGB1 NASP MCM6 GINS2 RAMP1 NCAPG TUBB4B C9orf16 GNG5 APOC1 HCST C1orf21 TYROBP CD247 FCER1G CD7 METTL1 CTSW CLIC3 CXXC5
24 CTLA4 TOX2 CD200 XAF1 RPS5 TOB1 HPGD CTSC NAMPT TNFRSF1B CFH IRAK3 JAML GAPDH C12orf75 EPSTI1 TSC22D3 PPP2R5C CTLA4 HIST1H1D MCM4 DEK CDCA7 ORC6 BRCA2 AREG RPL23A SAP18 RPL22 IKZF2 LITAF ADGRG1 AREG NFKBID GZMB GFOD1 CD7 PLSCR1 GSTP1
25 CCL4 CH25H BHLHE40-AS1 TMEM123 RPL3 AQP3 SELL DUSP4 LTA UGP2 NCR3 CA2 TNIP3 TNS3 CD74 OAS3 FOSB A2M-AS1 GAPDH TPX2 GINS2 H2AFZ CHEK1 MCM2 ECT2 CENPE RPL28 CIB1 PABPC1 FCRL6 PLAC8 CST7 ZEB2 ZFP36 SRGAP3 RGS16 CCL3 OAS1 IL2RB
26 FYB1 THADA PTPN13 DDX58 EEF1A1 RPL11 ICOS CD27 DNPH1 IFI44L FAM241A KIAA1211L TNFSF4 LAYN CD3G XAF1 ZFP36 GZMH HIST1H1B H2AFX CENPU HMGN2 MCM2 TTC38 CKAP2L UBE2S COX7C COPE HLA-DPB1 IL32 S100A4 HOPX AKR1C3 CEBPD CEBPD PARK7 LAT2 IFI35 LAT2
27 RUNX3 CD84 PDE7B CMPK2 RPL11 TNFAIP3 AC133644.2 ARID5B TIGIT SAMD9 CEBPD SLC4A10 CD63 NBL1 THEMIS USP18 PPP1R14B PDCD4 SMC4 SMC2 FEN1 SMC2 TYMS GZMB ARHGAP11B CKS1B CHCHD2 TRMT112 TPT1 NUCB2 SPON2 CD247 KLRB1 MAFF CD63 NDFIP2 GSTP1 KLRD1 CSF1
28 RNF19A BTLA NKG7 CD40LG RPS23 TNFRSF25 NABP1 CORO1B TFRC TNFRSF18 RPLP0 COL4A4 HLA-DQA1 MYO1E LINC02446 LY6E DUSP1 HLA-DPB1 SMC2 TMPO TUBB DNMT1 CDT1 PCNA ANLN TNFRSF18 RPL4 PSMB9 CD74 CD9 FCRL6 PRSS23 PTPN12 IGFBP2 NKG7 IFNG HOPX EPSTI1 SPRY1
29 PRF1 CORO1B IFNG SAMD9 RPL23 FXYD5 CD27 ENTPD1 F5 SAMD9L RPS12 NRP1 ID2 RGS2 HLA-DQB1 CMPK2 CD55 F2R BIRC5 CKS1B CHEK1 TPI1 MALAT1 PLEK ESPL1 HIST1H1B RPS24 BRK1 GZMH CD8A ADRB2 ABHD17A PRSS23 KDM6B TNFRSF18 COL6A3 CCL5 IRF7 CAPG
30 RGS1 METTL8 NR3C1 LY6E TPT1 RPL3 IKZF4 MAGEH1 PMAIP1 TBC1D4 ERN1 B3GALT5 FABP5 MYO7A LGALS1 NT5C3A CD69 SLAMF7 PKM UBE2S KNL1 ATAD2 HLA-DQA1 MCM5 NUF2 KPNA2 COX8A RBX1 IGHM KLRC4 ARL4C CTSW C1orf21 IER3 IFITM2 CXCL13 KLRF1 CEBPD ITGA1
31 GZMB CCDC50 LIMS1 GPR183 RPLP0 RPL9 PHTF2 S100A4 CCND2 BATF KIF5C PDZK1 LINC02446 AKAP5 HLA-DRB5 IFI44 GPR183 TNFSF9 RRM2 KNL1 SMC2 BRCA1 MYO7A HES6 KIF4A DLGAP5 RPS16 HSPA8 RPL8 MATK S1PR5 ZEB2 TXK CD160 CD300A TNFRSF18 SH2D1B XCL2 CLIC3
32 IL6R SMCO4 RNF19A USP18 RPL30 RPL32 PHACTR2 BIRC3 DUSP4 USP18 RPL13 EMID1 TNFRSF9 DGKH S100A6 HELZ2 TYROBP STK17A PMAIP1 TUBB4B MCM5 HNRNPAB MKI67 PRF1 KIFC1 CKS2 RPL35A S100A11 CD3E CRTAM KLF3 CHST2 S1PR5 REL GSTP1 RGCC PRF1 HERC5 LDLRAD4
33 CST7 CD4 RDH10 DDX60L RPS28 RPS3A RGS1 CD4 LINC02099 COX5A SPOCK2 HOXA5 TIGIT SNX9 FABP5 IRF7 FCER1G PKM TK1 HIST1H1C KIF11 SLBP HLA-DPA1 FEN1 BUB1 XCL2 RPL7 VAMP8 HLA-DQB1 S100A6 PATL2 PTPN12 CEBPD GNLY KRT81 CD82 ITGA1 NT5C3A KIR2DL4
34 SLC9A9 ARMH1 GAPDH TRIM22 LEF1 RPS6 DUSP16 MAF CCL20 SOX4 JAML LINGO4 ITGAE DUSP4 ANXA5 CD38 CCL3 COX5A JPT1 EZH2 TPX2 SNRPB SH3BGRL3 CDT1 KIF20B SMC2 RPS7 POMP WNK1 CAPG ANXA1 MYBL1 TRDC KLRF1 SPTSSB MALAT1 SRGAP3 IL2RB IGFBP2
35 BTLA PHACTR2 SNX9 HELZ2 RPS20 CD4 TRAC SPOCK2 SH2D2A IKZF2 TIGIT NECTIN2 LYST LINC01480 CCL4L2 TRIM22 CTSC EIF2AK2 ANP32B DEK HIST1H1B GZMB MCM3 ADGRG1 FBXO43 TMPO RPS11 EDF1 CD4 IFITM2 TGFBR3 TTC38 CTSW CLIC3 IGFBP2 RBPJ MCTP2 LY6E TXK
36 TP53INP1 TNFRSF25 IL6ST ANXA1 RPL38 RPL36 CARD16 IL1R1 ZBTB32 PMAIP1 RPL34 HOXA7 APOBEC3C ITGAE TMSB4X DDX58 RGCC ACP5 TMPO CDK1 GMNN COTL1 PTTG1 CCDC28B NCAPG2 CDK1 ATP5MG NDUFS5 HLA-DRA MAF CTSW S1PR5 GZMH TXK CAPN12 RANBP1 MAPK1 XAF1 LINC01871
37 TNFSF8 ARID5B CORO1B NT5C3A RPS16 RPL30 LINC01943 GLRX ID3 CTSC S100A6 PRAM1 CD27 RGS1 LGALS3 PARP14 AC020916.1 ADGRG1 CENPE BIRC5 MCM7 DNAJC9 MCM7 GNLY PRC1 CTSW MYL12B UQCR11 RPS13 CD4 S100A6 ITGB2 TTC38 MATK CXXC5 EIF5A CXXC5 XCL1 CD9
38 GZMH CTLA4 ITGA2 PLSCR1 RPS4X RPS8 S100A4 CARD16 SOX4 NT5C3A RPL32 B4GALNT1 HLA-DMA CD63 ITGA1 CD8A CTLA4 IFI35 UBE2C H2AFV RRM1 HLA-DRB1 RBPJ ASCL2 DIAPH3 KNL1 RPL39 LDHB CD3G GLUL TTC16 XBP1 ABHD17A KLRB1 IL18 DCTPP1 KLRC2 PRF1 SLC16A3
39 GNLY ZBED2 SMCO4 IFIH1 RPL36 RPS27A TNFRSF1B TYMP PTP4A3 TNFRSF9 HLA-DPB1 LPAR1 CRTAM SNAP47 TRGC2 DDX60 CD4 CSF1 CLSPN TK1 WDR76 CENPX CD8B PTPN12 SMC4 KRT86 UBB HNRNPA1 RPL29 PRMT9 A2M-AS1 CD160 PTGDR NKG7 ID2 ENO1 SAMD3 SAMD9 CCL5
40 METTL8 TNFRSF18 SAMSN1 IFI35 RPLP2 RPSA MAGEH1 MIR4435-2HG MIR155HG IFI44 TLE1 SLCO2A1 NKG7 CLNK APOBEC3C PPM1K DNAJB1 NDFIP2 CKS2 CENPM RRM2 IFNG PDCD1 PCLAF RPL10 PCLAF RPS21 ANP32B TNFRSF18 SCML4 TTC38 LAIR2 MYBL1 EGR2 PTGDR TPI1 LINC00996 MX2 NCR3
41 ADAM19 RILPL2 SRGN IRF7 RPS27 RPL13 CD4 DNPH1 ARID5B ARID5B TMIGD2 CD33 CTLA4 TNIP3 SIT1 IFIH1 MYBL1 LGALS3 ACTB ZWINT CENPF CHEK1 BRCA1 CDC45 RPS27 TYROBP RPL13 ABRACL LARP7 LYAR ASCL2 KLRB1 CMC1 IL2RB KLRB1 CAMK1 CD247 OASL PRMT9
42 CD2 AGFG1 ICOS DDX60 RPS2 EEF1B2 SPOCK2 PHACTR2 CD27 LINC01943 FKBP11 APOL4 RGS1 ADAMTS6 JAML RNF213 CD38 TMEM123 CALM3 HIST1H1E UBE2T CCL3 ITGA2 GMNN BRCA1 XCL1 ANP32B TRAPPC1 DUSP23 CTSA VCL TXK CHST2 MAP3K8 SH2D1B TNIP3 PTGDR USP18 CLNK
43 SPOCK2 CYSLTR1 KLRD1 OAS2 RPL12 RPL38 CCNG2 IL1R2 ENO1 GBP5 HLA-DRB1 RASSF8 ACP5 BCL2L11 LINC01871 DDX60L TNFRSF4 IL6ST ATAD2 FAM111B CDCA7 CDT1 CD63 CENPU RPL41 HMMR COX6C ATP5F1D DUSP2 IFITM3 LILRB1 PTGDR CCL4 IRF8 C1orf162 PIM3 NCAM1 EIF2AK2 CTSA
44 GZMA CD40LG CD109 PNPT1 SNHG8 RPLP1 FCMR PHLDA1 IKZF2 OAS3 RPS27A TLE1 SNAP47 NDFIP2 CD2 OAS2 TNFAIP3 CD247 ANP32E KIF11 CDKN3 RAN CCL5 C1orf21 RPL30 TMIGD2 RPS28 SNRPB SELL PTGDR FGR CD300A SH2D1B KRT81 PLAC8 NAMPT STARD3NL ERICH3 KLRB1
45 CD8B BHLHE40-AS1 CD3D PPM1K RPS27A RPL5 TLK1 STAM CD4 CD4 KLRD1 AL034397.3 ITM2A RAB27A CST7 PARP9 HLA-DRA TPI1 ICOS FABP5 SGO2 PKM CCL3 FGR RPL13 CDKN3 RPL36 SEPT7 KLF2 SPRY1 PLEKHG3 FGR CEP78 SRGN KIR2DL4 TNFRSF1B CAPN12 CCL3 CTSD
46 MATK IGFBP4 SLAMF1 PARP9 RPL19 RPS16 SMAP2 SLAMF1 LAIR2 RNF213 RUNX2 STAC SAMSN1 TNFRSF1B CCR5 PNPT1 JUN ANKRD28 CKS1B DLGAP5 USP1 ENO1 HLA-DMA LYN RPS27A H2AFX LDHB ZFAS1 CXCR4 CTLA4 COL6A2 GSAP XBP1 LAT2 MCTP2 ZNF593 TNFRSF18 IFI44 NDFIP2
47 CD84 ELMO1 RPL23A LGALS3BP PABPC1 CD52 LAYN LAIR2 SYNGR2 SPOCK2 RPL10 BCAS1 CCND2 ID2 CRTAM GZMK NCF1 CD226 PHF19 GTSE1 FANCI SNRPD1 GZMH KIR2DL3 RPS12 CD63 SUMO2 UQCRH DNAJC7 DDAH2 CES1 CTBP2 SAMD3 ZBTB16 ZNF683 HSPA5 AOAH CD7 ADGRG3
48 ARID5B MS4A6A AHI1 LTB RPS14 LMNA PIM2 F5 SERPINB9 ENTPD1 HLA-DPA1 CSF2 APOBEC3G LYST GRAP2 BST2 TRDC TMEM173 TNFRSF18 TMEM106C MT-ND4 NUDT1 APOBEC3G ITGAM PRR11 CCNB1 RPL14 SEC61G TRBC2 GUK1 LINC02384 GK5 MATK GRASP MMP25-AS1 FAM3C GZMA OAS3 CD247
49 CD63 ST8SIA1 ARID5B LAMP3 RPS21 RPL37 SAMHD1 TRAC NR4A1 CD27 RPS13 TOX2 CD74 AHI1 IDH2 SP110 HAVCR2 TPM4 UBE2S ANP32B ATAD5 PGAM1 ZNF683 MCM10 RPS14 GTSE1 RPS19 PPP1CA SERPINB6 CD28 COTL1 C1orf21 ITGB2 NR4A3 MAFF SLC7A5 APOBEC3G KLRC2 GZMB
50 HOPX MAGEH1 TYMP GZMA RPL35A RPS29 ICA1 PHTF2 VDR LY6E GZMB ALDOC DUSP4 ZEB2 GZMM C19orf66 MYADM USP18 RANBP1 HNRNPAB CENPE ATAD5 ANXA5 CD300A TPT1 BIRC5 ATP5F1E RHOA TYROBP CCR7 CD27 MTSS1 CD300A EGR3 XBP1 PGAM1 METRNL WARS CKLF
write_tsv(marker_sheet_extended, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_marker_sheet_full.tsv"))

write_tsv(marker_sheet_extended, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/supplementary_tables/", coi, "_marker_sheet_full.tsv"))

write_tsv(marker_tbl_extended_annotated, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_marker_table_annotated_full.tsv"))

write_tsv(marker_tbl_extended_annotated, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/supplementary_tables/", coi, "_marker_table_annotated_full.tsv"))
plot_data_sub_sub <- as_tibble(FetchData(seu_obj_sub_sub, c(myfeatures, "cluster_label_sub"))) %>% 
  left_join(meta_tbl, by = "sample") %>% 
  mutate(cell_id = colnames(seu_obj_sub_sub),
         cluster_label_sub = ordered(cluster_label_sub, levels = names(clrs$cluster_label_sub[[coi]])),
         ) %>%
  mutate(patient_id = str_remove_all(patient_id, "SPECTRUM-OV-"),
         tumor_supersite = ordered(tumor_supersite, levels = names(clrs$tumor_supersite)))
  
if (cell_sort == "CD45+") {
plot_data_sub_sub <- filter(plot_data_sub_sub, sort_parameters != "singlet, live, CD45-", !is.na(tumor_supersite))
}

if (cell_sort == "CD45-") {
plot_data_sub_sub <- filter(plot_data_sub_sub, sort_parameters != "singlet, live, CD45+", !is.na(tumor_supersite))
}

ggplot(plot_data_sub_sub, aes(umapharmony_1, umapharmony_2, color = cluster_label_sub)) + 
  common_layers_disc +
  # facet_wrap(~cluster_label_sub) +
  scale_color_manual(values = clrs$cluster_label_sub[[coi]]) +
  ggtitle("Sub cluster") 

4.1 add module scores and pathway scores

# signature_modules <- read_excel("_data/small/signatures/SPECTRUM v5 sub cluster markers.xlsx", sheet = 2, skip = 1, range = "M2:P100") %>%
#   gather(module, gene) %>%
#   na.omit() %>%
#   group_by(module) %>%
#   do(gene = c(.$gene)) %>%
#   {setNames(.$gene, .$module)}
# 
# signature_modules$ISG.module <- c("CCL5", "CXCL10", "IFNA1", "IFNB1", "ISG15", "IFI27L2", "SAMD9L")
# 
# ## compute expression module scores
# for (i in 1:length(signature_modules)) {
#   seu_obj_sub_sub <- AddModuleScore(seu_obj_sub_sub, features = signature_modules[i], name = names(signature_modules)[i])
#   seu_obj_sub_sub[[names(signature_modules)[i]]] <- seu_obj_sub_sub[[paste0(names(signature_modules)[i], "1")]]
#   seu_obj_sub_sub[[paste0(names(signature_modules)[i], "1")]] <- NULL
#   print(paste(names(signature_modules)[i], "DONE"))
# }
# 
# ## compute progeny scores
# progeny_list <- seu_obj_sub_sub@assays$RNA@data[VariableFeatures(seu_obj_sub_sub),] %>%
#   as.matrix %>%
#   progeny %>%
#   as.data.frame %>%
#   as.list
# 
# names(progeny_list) <- make.names(paste0(names(progeny_list), ".pathway"))
# 
# for (i in 1:length(progeny_list)) {
#   seu_obj_sub_sub <- AddMetaData(seu_obj_sub_sub,
#                                  metadata = progeny_list[[i]],
#                                  col.name = names(progeny_list)[i])
# }
# 
# write_rds(seu_obj_sub_sub, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered_annotated.rds"))

seu_obj_sub_sub <- read_rds(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered_annotated.rds"))

4.2 marker heatmap

marker_top_tbl <- marker_sheet_extended[,-1] %>% 
  slice(1:10) %>% 
  as.list %>% 
  .[!str_detect(names(.), "doublet|dissociated")] %>% 
  enframe("cluster_label_x", "gene") %>% 
  unnest(gene)

plot_data_markers <- as_tibble(FetchData(seu_obj_sub_sub, c("cluster_label_sub", myfeatures, unique(marker_top_tbl$gene)))) %>% 
  gather(gene, value, -c(1:(length(myfeatures)+1))) %>% 
  left_join(meta_tbl, by = "sample") %>% 
  mutate(patient_id = str_remove_all(patient_id, "SPECTRUM-OV-"),
         tumor_supersite = ordered(tumor_supersite, levels = names(clrs$tumor_supersite))) %>% 
  mutate(cluster_label_sub = ordered(cluster_label_sub, levels = names(clrs$cluster_label_sub[[coi]]))) %>% 
  group_by(cluster_label_sub, gene) %>% 
  summarise(value = mean(value, na.rm = T)) %>% 
  group_by(gene) %>% 
  mutate(value = scales::rescale(value)) %>% 
  left_join(marker_top_tbl, by = "gene") %>% 
  mutate(cluster_label_x = ordered(cluster_label_x, levels = rev(names(clrs$cluster_label_sub[[coi]]))))

ggplot(plot_data_markers) +
  geom_tile(aes(gene, cluster_label_sub, fill = value)) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  facet_wrap(~cluster_label_x, scales = "free") +
  scale_fill_gradientn(colors = viridis(9)) +
  labs(fill = "Scaled\nexpression") +
  theme(aspect.ratio = 1,
        axis.line = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank())

# ggsave(paste0("_fig/002_marker_heatmap_", coi, ".pdf"), width = nrow(marker_top_tbl)/6, height = 5)

4.3 composition

4.3.1 per site

comp_site_tbl <- plot_data_sub_sub %>%
  filter(!is.na(tumor_supersite)) %>% 
  group_by(cluster_label_sub, tumor_supersite) %>%
  tally %>%
  group_by(tumor_supersite) %>%
  mutate(nrel = n/sum(n)*100) %>%
  ungroup

pnrel_site <- ggplot(comp_site_tbl) +
  geom_bar(aes(tumor_supersite, nrel, fill = cluster_label_sub),
           stat = "identity", position = position_stack()) +
  scale_y_continuous(expand = c(0, 0)) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  labs(fill = "Cluster", y = "Fraction [%]", x = "") +
  scale_fill_manual(values = clrs$cluster_label_sub[[coi]])

pnabs_site <- ggplot(comp_site_tbl) +
  geom_bar(aes(tumor_supersite, n, fill = cluster_label_sub),
           stat = "identity", position = position_stack()) +
  scale_y_continuous(expand = c(0, 0)) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  labs(fill = "Cluster", y = "# cells", x = "") +
  scale_fill_manual(values = clrs$cluster_label_sub[[coi]])

plot_grid(pnabs_site, pnrel_site, ncol = 2, align = "h")

# ggsave(paste0("_fig/02_deep_dive_", coi, "_comp_site.pdf"), width = 8, height = 4)

4.3.2 per sample

comp_tbl_sample_sort <- plot_data_sub_sub %>% 
  group_by(tumor_subsite, tumor_supersite, patient_id, consensus_signature, therapy, cluster_label_sub) %>% 
  tally %>% 
  group_by(tumor_subsite, tumor_supersite, patient_id, consensus_signature, therapy) %>% 
  mutate(nrel = n/sum(n)*100,
         nsum = sum(n),
         log10n = log10(n),
         label_supersite = "Site",
         label_mutsig = "Signature",
         label_therapy = "Rx") %>% 
  ungroup %>% 
  arrange(desc(therapy), tumor_supersite) %>% 
  mutate(tumor_subsite_rx = paste0(tumor_subsite, "_", therapy)) %>% 
  mutate(tumor_subsite = ordered(tumor_subsite, levels = unique(tumor_subsite)),
         tumor_subsite_rx = ordered(tumor_subsite_rx, levels = unique(tumor_subsite_rx))) %>% 
  arrange(patient_id) %>% 
  mutate(label_patient_id = ifelse(as.logical(as.numeric(fct_inorder(as.character(patient_id)))%%2), "Patient1", "Patient2"))

sample_id_x_tbl <- plot_data_sub %>% 
  mutate(sort_short_x = cell_sort) %>% 
  distinct(patient_id, sort_short_x, tumor_subsite, therapy, sample) %>% 
  unite(sample_id_x, patient_id, sort_short_x, tumor_subsite, therapy) %>% 
  arrange(sample_id_x)

comp_tbl_sample_sort %>% 
  mutate(sort_short_x = cell_sort) %>% 
  unite(sample_id_x, patient_id, sort_short_x, tumor_subsite_rx) %>% 
  select(sample_id_x, cluster_label_sub, n, nrel, nsum) %>% 
  left_join(sample_id_x_tbl, by = "sample_id_x") %>% 
  write_tsv(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_subtype_compositions.tsv"))
ybreaks <- c(500, 1000, 2000, 4000, 6000, 8000, 10000, 15000, 20000)

max_cells_per_sample <- max(comp_tbl_sample_sort$nsum)
ymaxn <- ybreaks[max_cells_per_sample < ybreaks][1]

comp_plot_wrapper <- function(y = "nrel", switch = NULL) {
  if (y == "nrel") ylab <- paste0("Fraction\nof cells [%]")
  if (y == "n") ylab <- paste0("Number\nof cells")
  p <- ggplot(comp_tbl_sample_sort, 
              aes_string("tumor_subsite_rx", y, fill = "cluster_label_sub")) + 
    facet_grid(~patient_id, space = "free", scales = "free", switch = switch) +
    coord_cartesian(clip = "off") + 
    scale_fill_manual(values = clrs$cluster_label_sub[[coi]]) + 
    theme(axis.text.x = element_blank(),
          axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5, 
                                      margin = margin(0, -0.4, 0, 0, unit = "npc")),
          axis.ticks.x = element_blank(),
          axis.title.x = element_blank(),
          axis.line.x = element_blank(),
          strip.text.y = element_blank(),
          strip.text.x = element_blank(),
          strip.background.y = element_blank(),
          strip.background.x = element_blank(),
          plot.margin = margin(0, 0, 0, 0, "npc")) + 
    labs(x = "",
         y = ylab) +
    guides(fill = FALSE)
  if (y == "nrel") p <- p + 
    geom_bar(stat = "identity") +
    scale_y_continuous(expand = c(0, 0), 
                       breaks = c(0, 50, 100), 
                       labels = c("0", "50", "100"))
  if (y == "n") p <- p + 
    geom_bar(stat = "identity", position = position_stack()) +
    scale_y_continuous(expand = c(0, 0), 
                       breaks = c(0, ymaxn/2, ymaxn),
                       limits = c(0, ymaxn),
                       labels = c("", ymaxn/2, ymaxn)) +
    expand_limits(y = c(0, ymaxn)) +
    theme(panel.grid.major.y = element_line(linetype = 1, color = "grey90", size = 0.5))
  return(p)
} 

common_label_layers <- list(
  geom_tile(color = "white", size = 0.15),
  theme(axis.text.x = element_blank(),
        axis.ticks = element_blank(),
        axis.title.x = element_blank(),
        axis.line.x = element_blank(),
        strip.text = element_blank(),
        strip.background = element_blank(),
        plot.margin = margin(0, 0, 0, 0, "npc")),
  scale_y_discrete(expand = c(0, 0)),
  labs(y = ""),
  guides(fill = FALSE),
  facet_grid(~patient_id, 
             space = "free", scales = "free")
)

comp_label_site <- ggplot(distinct(comp_tbl_sample_sort, tumor_subsite_rx, therapy, tumor_supersite, label_supersite, patient_id), 
       aes(tumor_subsite_rx, label_supersite, 
           fill = tumor_supersite)) + 
  scale_fill_manual(values = clrs$tumor_supersite) +
  common_label_layers

comp_label_rx <- ggplot(distinct(comp_tbl_sample_sort, tumor_subsite_rx, therapy, tumor_supersite, label_therapy, consensus_signature, patient_id), 
       aes(tumor_subsite_rx, label_therapy, 
           fill = therapy)) + 
  scale_fill_manual(values = c(`post-Rx` = "gold3", `pre-Rx` = "steelblue")) +
  common_label_layers

comp_label_mutsig <- ggplot(distinct(comp_tbl_sample_sort, tumor_subsite_rx, therapy, tumor_supersite, label_mutsig, consensus_signature, patient_id), 
       aes(tumor_subsite_rx, label_mutsig, 
           fill = consensus_signature)) + 
  scale_fill_manual(values = clrs$consensus_signature) +
  common_label_layers

patient_label_tbl <- distinct(comp_tbl_sample_sort, patient_id, .keep_all = T)

comp_label_patient_id <- ggplot(patient_label_tbl, aes(tumor_subsite_rx, label_patient_id)) + 
  scale_fill_manual(values = clrs$consensus_signature) +
  geom_text(aes(tumor_subsite_rx, label_patient_id, label = patient_id)) +
  facet_grid(~patient_id, 
             space = "free", scales = "free") +
  coord_cartesian(clip = "off") + 
  theme_void() +
  theme(strip.text = element_blank(),
        strip.background = element_blank(),
        plot.margin = margin(0, 0, 0, 0, "npc"))

hist_plot_wrapper <- function(x = "nrel") {
  if (x == "nrel") {
    xlab <- paste0("Fraction of cells [%]")
    bw <- 5
  }
  if (x == "log10n") {
    xlab <- paste0("Number of cells")
    bw <- 0.2
  }
  p <- ggplot(comp_tbl_sample_sort) +
    ggridges::geom_density_ridges(
      aes_string(x, "cluster_label_sub", fill = "cluster_label_sub"), color = "black",
      stat = "binline", binwidth = bw, scale = 3) +
    facet_grid(label_supersite~., 
               space = "free", scales = "free") +
    scale_fill_manual(values = clrs$cluster_label_sub[[coi]]) + 
    theme(axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          axis.title.y = element_blank(),
          axis.line.y = element_blank(),
          strip.text = element_blank(),
          strip.background = element_blank(),
          plot.margin = margin(0, 0, 0, 0, "npc")) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_discrete(expand = c(0, 0)) +
    guides(fill = FALSE) +
    labs(x = xlab)
  if (x == "log10n") p <- p + expand_limits(x = c(0, 3)) + 
    scale_x_continuous(expand = c(0, 0), 
                       labels = function(x) parse(text = paste("10^", x)))
  return(p)
}

pcomp1 <- comp_plot_wrapper("n")
pcomp2 <- comp_plot_wrapper("nrel")

pcomp_grid <- plot_grid(comp_label_patient_id, 
                        pcomp1, pcomp2, 
                        comp_label_site, comp_label_rx, comp_label_mutsig,
                        ncol = 1, align = "v", 
                        rel_heights = c(0.15, 0.33, 0.33, 0.06, 0.06, 0.06))

phist1 <- hist_plot_wrapper("log10n")

pcomp_hist_grid <- ggdraw() +
  draw_plot(pcomp_grid, x = 0.01, y = 0, width = 0.85, height = 1) +
  draw_plot(phist1, x = 0.87, y = 0.05, width = 0.12, height = 0.8)

pcomp_hist_grid

# ggsave(paste0("_fig/02_composition_v6_",coi,".pdf"), pcomp_hist_grid, width = 10, height = 2)

4.3.3 site specific cluster enrichment

comp_tbl_z <- comp_tbl_sample_sort %>% 
  filter(therapy == "pre-Rx",
         !(tumor_supersite %in% c("Ascites", "Other"))) %>% 
  group_by(patient_id, cluster_label_sub) %>% 
  arrange(patient_id, cluster_label_sub, nrel) %>% 
  mutate(rank = row_number(nrel),
         z_rank = scales::rescale(rank)) %>% 
  mutate(mean_nrel = mean(nrel, na.rm = T),
         sd_nrel = sd(nrel, na.rm = T),
         z_nrel = (nrel - mean_nrel) / sd_nrel) %>% 
  ungroup()

ggplot(comp_tbl_z) +
  geom_boxplot(aes(tumor_supersite, z_nrel, color = tumor_supersite), 
               outlier.shape = NA) +
  geom_point(aes(tumor_supersite, z_nrel, color = tumor_supersite), position = position_jitter(), size = 0.1) +
  facet_grid(~cluster_label_sub, scales = "free_x", space = "free_x") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        aspect.ratio = 5) +
  scale_color_manual(values = clrs$tumor_supersite)

ggplot(comp_tbl_z) +
  geom_boxplot(aes(tumor_supersite, z_rank, color = tumor_supersite), 
               outlier.shape = NA) +
  geom_point(aes(tumor_supersite, z_rank, color = tumor_supersite), position = position_jitter(), size = 0.1) +
  facet_grid(~cluster_label_sub, scales = "free_x", space = "free_x") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        aspect.ratio = 5) +
  scale_color_manual(values = clrs$tumor_supersite)

## sub cluster CD8 cells
# seu_obj_cd8 <- seu_obj_sub_sub %>%
#   subset(subset = cluster_label_sub %in% grep("^CD8.T", unique(seu_obj_sub_sub$cluster_label_sub), value = T)) %>%
#   RunUMAP(dims = 1:50, reduction = "harmony", reduction.name = "umapharmony", reduction.key = "umapharmony_")
# 
# write_rds(seu_obj_cd8, "/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T_processed.rds")

seu_obj_cd8 <- read_rds("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T_processed.rds")
root_cell <- "SPECTRUM-OV-041_S1_CD45P_RIGHT_FALLOPIAN_TUBE_ACCCTCATCCAAGCAT"
seu_obj_cd8_sub <- subset(seu_obj_cd8, cells = c(root_cell, colnames(seu_obj_cd8)[colnames(seu_obj_cd8)!=root_cell][-1]))

dc_obj <- DiffusionMap(seu_obj_cd8_sub@reductions$harmony@cell.embeddings, k = 100)
dc_mat <- dc_obj@eigenvectors
colnames(dc_mat) <- paste0("DC_", 1:ncol(dc_mat))
seu_obj_cd8_sub[["DC"]] <- CreateDimReducObject(embeddings = dc_mat, key = "DC_", assay = DefaultAssay(seu_obj_cd8_sub))

write_rds(seu_obj_cd8_sub, "/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T_processed_filtered.rds")

dpt_obj <- DPT(dc_obj)
write_rds(dpt_obj, "/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T_sub_dpt.rds")
seu_obj_cd8_sub$DPT1 <- dpt_obj$DPT1

write_rds(seu_obj_cd8_sub, "/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T_processed_filtered.rds")
seu_obj_cd8_sub <- read_rds("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T_processed_filtered.rds")
# plot_data_sub_sub <- as_tibble(FetchData(seu_obj_sub_sub, c(myfeatures, "cluster_label"))) %>% 
#   left_join(meta_tbl, by = "sample") %>% 
#   mutate(patient_id = str_remove_all(patient_id, "SPECTRUM-OV-"),
#          tumor_supersite = ordered(tumor_supersite, levels = names(clrs$tumor_supersite))) %>% 
#   mutate(cell_id = colnames(seu_obj_sub_sub),
#          cluster_label = ordered(cluster_label, levels = my_subtypes),
#          )
#   
# if (cell_sort == "CD45+") {
# plot_data_sub_sub <- filter(plot_data_sub_sub, sort_parameters != "singlet, live, CD45-", !is.na(tumor_supersite))
# }
# 
# if (cell_sort == "CD45-") {
# plot_data_sub_sub <- filter(plot_data_sub_sub, sort_parameters != "singlet, live, CD45+", !is.na(tumor_supersite))
# }
# 
# ggplot(plot_data_sub_sub, aes(umapharmony_1, umapharmony_2, color = cluster_label)) + 
#   common_layers_disc +
#   ggtitle("Sub cluster") +
#   #facet_wrap(~cluster_label) +
#   scale_color_manual(values = clrs$cluster_label[[coi]])
# 
# ggplot(plot_data_sub_sub, aes(umapharmony_1, umapharmony_2, color = patient_id_short)) + 
#   common_layers_disc +
#   ggtitle("Patient") +
#   # facet_wrap(~therapy) +
#   scale_color_manual(values = clrs$patient_id_short)
# 
# ggplot(plot_data_sub_sub, aes(umapharmony_1, umapharmony_2, color = tumor_supersite)) + 
#   # geom_point(aes(umapharmony_1, umapharmony_2), 
#   #            color = "grey90", size = 0.01, 
#   #            data = select(plot_data_sub_sub, -tumor_supersite)) +
#   common_layers_disc +
#   ggtitle("Site") +
#   # facet_wrap(~therapy) +
#   scale_color_manual(values = clrs$tumor_supersite)
# 
# write_tsv(select(plot_data_sub_sub, cell_id, everything(), -umapharmony_1, -umapharmony_2, -contains("RNA_")), paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_embedding_sub.tsv"))

4.3.4 UMAP & DC

plot_data_cd8_sub <- as_tibble(FetchData(seu_obj_cd8_sub, c(myfeatures, "cluster_label_sub", "DC_1", "DC_2"))) %>% 
  left_join(meta_tbl, by = "sample") %>% 
  mutate(patient_id = str_remove_all(patient_id, "SPECTRUM-OV-"),
         tumor_supersite = ordered(tumor_supersite, levels = names(clrs$tumor_supersite))) %>% 
  mutate(cell_id = colnames(seu_obj_cd8_sub),
         cluster_label_sub = ordered(cluster_label_sub, levels = names(clrs$cluster_label_sub[[coi]])))
  
if (cell_sort == "CD45+") {
plot_data_cd8_sub <- filter(plot_data_cd8_sub, sort_parameters != "singlet, live, CD45-", !is.na(tumor_supersite))
}

if (cell_sort == "CD45-") {
plot_data_cd8_sub <- filter(plot_data_cd8_sub, sort_parameters != "singlet, live, CD45+", !is.na(tumor_supersite))
}

ggplot(plot_data_cd8_sub, aes(umapharmony_1, umapharmony_2, color = cluster_label_sub)) + 
  common_layers_disc +
  ggtitle("Sub cluster") +
  #facet_wrap(~cluster_label_sub) +
  scale_color_manual(values = clrs$cluster_label_sub[[coi]])

ggplot(plot_data_cd8_sub, aes(umapharmony_1, umapharmony_2, color = patient_id_short)) + 
  common_layers_disc +
  ggtitle("Patient") +
  # facet_wrap(~therapy) +
  scale_color_manual(values = clrs$patient_id_short)

ggplot(plot_data_cd8_sub, aes(umapharmony_1, umapharmony_2, color = tumor_supersite)) + 
  common_layers_disc +
  ggtitle("Site") +
  # facet_wrap(~therapy) +
  scale_color_manual(values = clrs$tumor_supersite)

ggplot(plot_data_cd8_sub, aes(DC_1, DC_2, color = cluster_label_sub)) + 
  common_layers_disc +
  ggtitle("Sub cluster") +
  #facet_wrap(~cluster_label_sub) +
  scale_color_manual(values = clrs$cluster_label_sub[[coi]])

ggplot(plot_data_cd8_sub, aes(DC_1, DC_2, color = patient_id_short)) + 
  common_layers_disc +
  ggtitle("Patient") +
  # facet_wrap(~therapy) +
  scale_color_manual(values = clrs$patient_id_short)

ggplot(plot_data_cd8_sub, aes(DC_1, DC_2, color = tumor_supersite)) + 
  common_layers_disc +
  ggtitle("Site") +
  # facet_wrap(~therapy) +
  scale_color_manual(values = clrs$tumor_supersite)

5 session info

devtools::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 3.6.2 (2019-12-12)
##  os       Debian GNU/Linux 10 (buster)
##  system   x86_64, linux-gnu           
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  ctype    en_US.UTF-8                 
##  tz       Etc/UTC                     
##  date     2021-02-07                  
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package              * version    date       lib
##  abind                  1.4-5      2016-07-21 [2]
##  ape                    5.3        2019-03-17 [2]
##  assertthat             0.2.1      2019-03-21 [2]
##  backports              1.1.10     2020-09-15 [1]
##  bibtex                 0.4.2.2    2020-01-02 [2]
##  Biobase                2.46.0     2019-10-29 [2]
##  BiocGenerics           0.32.0     2019-10-29 [2]
##  BiocParallel           1.20.1     2019-12-21 [2]
##  bitops                 1.0-6      2013-08-17 [2]
##  boot                   1.3-24     2019-12-20 [3]
##  broom                  0.7.2      2020-10-20 [1]
##  callr                  3.4.2      2020-02-12 [1]
##  car                    3.0-8      2020-05-21 [1]
##  carData                3.0-4      2020-05-22 [1]
##  caTools                1.17.1.4   2020-01-13 [2]
##  cellranger             1.1.0      2016-07-27 [2]
##  class                  7.3-15     2019-01-01 [3]
##  cli                    2.0.2      2020-02-28 [1]
##  cluster                2.1.0      2019-06-19 [3]
##  codetools              0.2-16     2018-12-24 [3]
##  colorblindr          * 0.1.0      2020-01-13 [2]
##  colorspace           * 1.4-2      2019-12-29 [2]
##  cowplot              * 1.0.0      2019-07-11 [2]
##  crayon                 1.3.4      2017-09-16 [1]
##  curl                   4.3        2019-12-02 [2]
##  data.table             1.12.8     2019-12-09 [2]
##  DBI                    1.1.0      2019-12-15 [2]
##  dbplyr                 2.0.0      2020-11-03 [1]
##  DelayedArray           0.12.2     2020-01-06 [2]
##  DEoptimR               1.0-8      2016-11-19 [1]
##  desc                   1.2.0      2018-05-01 [2]
##  destiny              * 3.0.1      2020-01-16 [1]
##  devtools               2.2.1      2019-09-24 [2]
##  digest                 0.6.25     2020-02-23 [1]
##  dplyr                * 1.0.2      2020-08-18 [1]
##  e1071                  1.7-3      2019-11-26 [1]
##  ellipsis               0.3.1      2020-05-15 [1]
##  evaluate               0.14       2019-05-28 [2]
##  fansi                  0.4.1      2020-01-08 [2]
##  farver                 2.0.3      2020-01-16 [1]
##  fitdistrplus           1.0-14     2019-01-23 [2]
##  forcats              * 0.5.0      2020-03-01 [1]
##  foreign                0.8-74     2019-12-26 [3]
##  formattable            0.2.0.1    2016-08-05 [1]
##  fs                     1.5.0      2020-07-31 [1]
##  future                 1.15.1     2019-11-25 [2]
##  future.apply           1.4.0      2020-01-07 [2]
##  gbRd                   0.4-11     2012-10-01 [2]
##  gdata                  2.18.0     2017-06-06 [2]
##  generics               0.0.2      2018-11-29 [2]
##  GenomeInfoDb           1.22.0     2019-10-29 [2]
##  GenomeInfoDbData       1.2.2      2020-01-14 [2]
##  GenomicRanges          1.38.0     2019-10-29 [2]
##  ggplot.multistats      1.0.0      2019-10-28 [1]
##  ggplot2              * 3.3.2      2020-06-19 [1]
##  ggrepel                0.8.1      2019-05-07 [2]
##  ggridges               0.5.2      2020-01-12 [2]
##  ggthemes               4.2.0      2019-05-13 [1]
##  globals                0.12.5     2019-12-07 [2]
##  glue                   1.3.2      2020-03-12 [1]
##  gplots                 3.0.1.2    2020-01-11 [2]
##  gridExtra              2.3        2017-09-09 [2]
##  gtable                 0.3.0      2019-03-25 [2]
##  gtools                 3.8.1      2018-06-26 [2]
##  haven                  2.3.1      2020-06-01 [1]
##  hexbin                 1.28.0     2019-11-11 [2]
##  hms                    0.5.3      2020-01-08 [1]
##  htmltools              0.5.1.1    2021-01-22 [1]
##  htmlwidgets            1.5.1      2019-10-08 [2]
##  httr                   1.4.2      2020-07-20 [1]
##  ica                    1.0-2      2018-05-24 [2]
##  igraph                 1.2.5      2020-03-19 [1]
##  IRanges                2.20.2     2020-01-13 [2]
##  irlba                  2.3.3      2019-02-05 [2]
##  jsonlite               1.7.1      2020-09-07 [1]
##  KernSmooth             2.23-16    2019-10-15 [3]
##  knitr                  1.26       2019-11-12 [2]
##  knn.covertree          1.0        2019-10-28 [1]
##  labeling               0.3        2014-08-23 [2]
##  laeken                 0.5.1      2020-02-05 [1]
##  lattice                0.20-38    2018-11-04 [3]
##  lazyeval               0.2.2      2019-03-15 [2]
##  leiden                 0.3.1      2019-07-23 [2]
##  lifecycle              0.2.0      2020-03-06 [1]
##  listenv                0.8.0      2019-12-05 [2]
##  lmtest                 0.9-37     2019-04-30 [2]
##  lsei                   1.2-0      2017-10-23 [2]
##  lubridate              1.7.9.2    2020-11-13 [1]
##  magrittr             * 2.0.1      2020-11-17 [1]
##  MASS                   7.3-51.5   2019-12-20 [3]
##  Matrix                 1.2-18     2019-11-27 [3]
##  matrixStats            0.56.0     2020-03-13 [1]
##  memoise                1.1.0      2017-04-21 [2]
##  metap                  1.2        2019-12-08 [2]
##  mnormt                 1.5-5      2016-10-15 [2]
##  modelr                 0.1.8      2020-05-19 [1]
##  multcomp               1.4-12     2020-01-10 [2]
##  multtest               2.42.0     2019-10-29 [2]
##  munsell                0.5.0      2018-06-12 [2]
##  mutoss                 0.1-12     2017-12-04 [2]
##  mvtnorm                1.0-12     2020-01-09 [2]
##  nlme                   3.1-143    2019-12-10 [3]
##  nnet                   7.3-12     2016-02-02 [3]
##  npsurv                 0.4-0      2017-10-14 [2]
##  numDeriv               2016.8-1.1 2019-06-06 [2]
##  openxlsx               4.1.5      2020-05-06 [1]
##  pbapply                1.4-2      2019-08-31 [2]
##  pcaMethods             1.78.0     2019-10-29 [2]
##  pillar                 1.4.6      2020-07-10 [1]
##  pkgbuild               1.0.6      2019-10-09 [2]
##  pkgconfig              2.0.3      2019-09-22 [1]
##  pkgload                1.0.2      2018-10-29 [2]
##  plotly                 4.9.1      2019-11-07 [2]
##  plotrix                3.7-7      2019-12-05 [2]
##  plyr                   1.8.5      2019-12-10 [2]
##  png                    0.1-7      2013-12-03 [2]
##  prettyunits            1.1.1      2020-01-24 [1]
##  processx               3.4.2      2020-02-09 [1]
##  progeny              * 1.11.3     2020-10-22 [1]
##  proxy                  0.4-24     2020-04-25 [1]
##  ps                     1.3.2      2020-02-13 [1]
##  purrr                * 0.3.4      2020-04-17 [1]
##  R.methodsS3            1.7.1      2016-02-16 [2]
##  R.oo                   1.23.0     2019-11-03 [2]
##  R.utils                2.9.2      2019-12-08 [2]
##  R6                     2.4.1      2019-11-12 [1]
##  ranger                 0.12.1     2020-01-10 [1]
##  RANN                   2.6.1      2019-01-08 [2]
##  rappdirs               0.3.1      2016-03-28 [2]
##  RColorBrewer           1.1-2      2014-12-07 [2]
##  Rcpp                   1.0.4      2020-03-17 [1]
##  RcppAnnoy              0.0.16     2020-03-08 [1]
##  RcppEigen              0.3.3.7.0  2019-11-16 [2]
##  RcppHNSW               0.2.0      2019-09-20 [2]
##  RcppParallel           4.4.4      2019-09-27 [2]
##  RCurl                  1.98-1.1   2020-01-19 [1]
##  Rdpack                 0.11-1     2019-12-14 [2]
##  readr                * 1.4.0      2020-10-05 [1]
##  readxl               * 1.3.1      2019-03-13 [2]
##  remotes                2.1.0      2019-06-24 [2]
##  reprex                 0.3.0      2019-05-16 [2]
##  reshape2               1.4.3      2017-12-11 [2]
##  reticulate             1.14       2019-12-17 [2]
##  rio                    0.5.16     2018-11-26 [1]
##  rlang                  0.4.8      2020-10-08 [1]
##  rmarkdown              2.0        2019-12-12 [2]
##  robustbase             0.93-6     2020-03-23 [1]
##  ROCR                   1.0-7      2015-03-26 [2]
##  rprojroot              1.3-2      2018-01-03 [2]
##  RSpectra               0.16-0     2019-12-01 [2]
##  rstudioapi             0.11       2020-02-07 [1]
##  rsvd                   1.0.3      2020-02-17 [1]
##  Rtsne                  0.15       2018-11-10 [2]
##  rvest                  0.3.6      2020-07-25 [1]
##  S4Vectors              0.24.2     2020-01-13 [2]
##  sandwich               2.5-1      2019-04-06 [2]
##  scales                 1.1.0      2019-11-18 [2]
##  scatterplot3d          0.3-41     2018-03-14 [1]
##  sctransform            0.2.1      2019-12-17 [2]
##  SDMTools               1.1-221.2  2019-11-30 [2]
##  sessioninfo            1.1.1      2018-11-05 [2]
##  Seurat               * 3.1.2      2019-12-12 [2]
##  SingleCellExperiment   1.8.0      2019-10-29 [2]
##  smoother               1.1        2015-04-16 [1]
##  sn                     1.5-4      2019-05-14 [2]
##  sp                     1.4-2      2020-05-20 [1]
##  stringi                1.5.3      2020-09-09 [1]
##  stringr              * 1.4.0      2019-02-10 [1]
##  SummarizedExperiment   1.16.1     2019-12-19 [2]
##  survival               3.1-8      2019-12-03 [3]
##  testthat               2.3.2      2020-03-02 [1]
##  TFisher                0.2.0      2018-03-21 [2]
##  TH.data                1.0-10     2019-01-21 [2]
##  tibble               * 3.0.4      2020-10-12 [1]
##  tidyr                * 1.1.2      2020-08-27 [1]
##  tidyselect             1.1.0      2020-05-11 [1]
##  tidyverse            * 1.3.0      2019-11-21 [2]
##  tsne                   0.1-3      2016-07-15 [2]
##  TTR                    0.23-6     2019-12-15 [1]
##  usethis                1.5.1      2019-07-04 [2]
##  uwot                   0.1.5      2019-12-04 [2]
##  vcd                    1.4-7      2020-04-02 [1]
##  vctrs                  0.3.5      2020-11-17 [1]
##  VIM                    6.0.0      2020-05-08 [1]
##  viridis              * 0.5.1      2018-03-29 [2]
##  viridisLite          * 0.3.0      2018-02-01 [2]
##  withr                  2.3.0      2020-09-22 [1]
##  xfun                   0.12       2020-01-13 [2]
##  xml2                   1.3.2      2020-04-23 [1]
##  xts                    0.12-0     2020-01-19 [1]
##  XVector                0.26.0     2019-10-29 [2]
##  yaml                   2.2.1      2020-02-01 [1]
##  zip                    2.0.4      2019-09-01 [1]
##  zlibbioc               1.32.0     2019-10-29 [2]
##  zoo                    1.8-7      2020-01-10 [2]
##  source                                 
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  Bioconductor                           
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Github (clauswilke/colorblindr@1ac3d4d)
##  R-Forge (R 3.6.2)                      
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  Bioconductor                           
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Github (saezlab/progeny@94bea1c)       
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.3)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
## 
## [1] /home/uhlitzf/R/lib
## [2] /usr/local/lib/R/site-library
## [3] /usr/local/lib/R/library